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Abstract 

We develop a model in the framework of nuclear fragmentation 
at thermodynamic equilibrium which can be mapped onto an Ising 
model with constant magnetization. We work out the thermodynamic 
properties of the model as well as the properties of the fragment size 
distributions. We show that two types of phase transitions can be 
found for high density systems. They merge into a unique transition 
at low density. An analysis of the critical exponents which characterize 
observables for different densities in the thermodynamic limit shows 
that these transitions look like continuous second order transitions. 



PACS: 25.70.Pq, 64.60.Cn, 64.60.Fr. 
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1 Introduction 



The fragmentation of excited nuclei is a complex, time-dependent process. 
There exists however experimental evidence if no proof that the decaying 
system can be treated as a thermodynamically equilibrated assembly of in- 
teracting particles, at least in specific experimental circumstances [Q. If this 
is the case, the fragmentation process evolves in time by following trajecto- 
ries in the (p, T) phase diagram (p is the density and T the temperature). 
It starts at some initial point and ends at a final point when the freeze-out 
stage is reached. During the time interval in which this goes on, the initial 
excited system breaks into fragments and (or) evaporates particles. 

As far as such a scenario corresponds to a physical realization of the pro- 
cess, one is left with two types of experimentally relevant physical observ- 
ables. The first set concerns the thermodynamic properties of the expanding 
system (energy, temperature, pressure and correlations between these quan- 
tities, possible existence of phase transitions). The second set concerns final 
fragment size distributions and related observables. 

The latter observables have been rather extensively studied in numer- 
ous experiments, showing remarkable properties |@] which were analyzed by 
means of different types of models. These analyses lead to the conclusion 
that the system may undergo a second order transition which lies in the 
universality class of bond percolation in three dimensions [^. 

The thermodynamic properties of fragmenting nuclei are less well known 
from the experimental point of view. There are indications for the existence 
of two phases which could be analogous to liquid and gas in macroscopic 
systems. The order of the transition between the two phases is not clearly 
established either |, ^. Successful multifragmentation models []^, |^, |^ 
which are able to reproduce a large amount of experimental data predict a 
first order transition. 

Running parallel to the experimental efforts in order to characterize the 
behaviour of the system, realistic theoretical models should aim and be able 
to describe both aspects, the thermodynamic properties as well as the char- 
acteristics of the fragment size distributions in a consistent way. This has 
already been attempted in recent work in the framework of the lattice gas 
model (LGM) |lTl, |l^ and to some extents by means of homogeneous 
cellular models . 



It is the aim of the present paper to present a detailed analysis of both 
properties by means of a simple though quite realistic model, to identify the 
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phases in which the system can exist and the order of the corresponding 
phase transitions. 

The paper is arranged in the following way. In section 2 we describe the 
model. We emphasize its kinship with other models and its specific features. 
In section 3 we first define the observables which we consider in the sequel 
and work out the phase diagram which links the density to the temperature. 
We then analyze the properties of the phase transitions which appear in the 
system in section 4. We summarize, conclude, present and discuss further 
possible developments in section 5. 

2 Spin models and fragmentation 

2.1 The lattice gas model and the Ising model 

The lattice gas model (LGM) describes the motion of N = p ■ V particles 
inside a three-dimensional cubic lattice, where p is the average density, < 
p < 1, and V the number of sites. Particles which occupy sites in this 
volume interact only when they are nearest-neighbours. If = 0, 1 are the 
site occupation numbers (empty/occupied), and p is the chemical potential, 
the Hamiltonian reads 

Hv = -4:J ^ niUj + p'^nj, (1) 
{i,j)cv j&v 

where {ij) are nearest neighbour cells. 

The formulation presented in eq. (|l]) corresponds to the grand canonical 
ensemble. The microcanonical ensemble is the most natural ensemble to be 
used when describing nuclear fragmentation since nuclei are closed systems 
with fixed energy and number of particles. It is usually assumed that the 
thermodynamic results are the same independently of the formalism used. 
Of course this is not true with finite systems, because microcanonical, canon- 
ical and grand canonical ensembles differ on the level of fluctuations. But it 
has recently been argued [|15| that even in the thermodynamic limit, micro- 
canonical and (grand) canonical ensembles are not equivalent (for example, at 
a first order transition). We shall explicitely show in this paper that in fact 
the canonical and grand canonical formulations of the LGM are not equiva- 
lent. This happens because there are points in the phase diagram which are 
accessible when using one formalism, but are not defined (they are not points 



3 



of thermodynamic equilibrium) when using another one. As we will see, this 
causes the order of the transition to be different in the two formalisms. This 
fact leads us naturally to ask about the behaviour of the system in the mi- 
crocanonical ensemble. This is an open problem, but as a first step we will 
consider here only the detailed study of the canonical formulation to com- 
pare it with the well-known grand canonical version of this model. In fact, 
for sytems with finite volume, the microcanonical ensemble turns out to be 



very efficient in the study of phase transitions [jT6|, [17 . 

This simple model can be numerically studied, its grand partition function 
being proportional to the canonical partition function of the Ising model with 



an external magnetic field [jT8|. In fact, setting 

a, = 2n, - 1, (2) 

the Hamiltonian (|l]) becomes 

Ry = -J ^ ^.QT^. + (^^_6jj^CT^. + const (r) (3) 

{i,j)<zv jev 

= —J cTjaj — h (Tj + const (l^), 

{i,j)cv jdV 

where we reinterpret = ±1 as a spin variable, and identify a spin up with 
an occupied site, and a spin down with an empty site. We introduce 

M^^Y.^, = ^{2N-V) = 2p-l. (4) 

j 

M is called the magnetization in the sequel. As it can be seen, the initial 
model defined through ([1|) can indeed be mapped on an Ising model with a 
linear term corresponding to a fixed magnetic field h. 

It is well known that the Ising model with magnetic field h presents a 
discontinuity in the magnetization M{h,T) for T < Tc ai h = 0, where Tc 
is the critical temperature. That is, in the plane {h, T) there exists a line of 
first-order transitions {h = 0, < T < Tc), and a second order transition at 
{h = 0,Tc). Let us see how this is refiected in the case of the LGM. Since 
h plays essentially the role of the chemical potential in the grand canonical 
formulation of the LGM, there is a first order transition in the phase diagram 
(p, T) of the LGM, which separates the homogeneous phase from the zone 
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with the two "hquid" and "gas" phases. On this hne, there exists a critical 
point at p = 0.5, T = Tc. 

As described below, we aim to introduce a similar type of model in order 
to describe a finite, excited system of classical particles which interact by 
means of a short range potential. We characterize the internal structure of 
the system by means of a proper definition of physical droplets. Droplets are 
clusters of particles in which each particle has at least one nearest neighbour. 
They exist because of two physical reasons: the attractive nearest-neighbour 
interaction and the density of particles, which leads to usual percolation when 
this density is large enough. However, such clusters are not always stable, 
because they can break at zones of small connectivity owing to the kinetic 
energy of the particles. The simple model (equivalently, the LGM) does not 
take the kinetic energy into account. We can introduce this energy by means 
of another definition which takes the energetics of the fragmentation into 
account [|1^, |ll|. Let us consider the stability of a cluster of A particles. If 
the loss of a particle implies the break of a single bond, then the excess of 
energy will be 

where m is the mass of the particles of the LGM. The cluster is a stable 
droplet if this excess is positive: 

^ - J < 0. (6) 

This suggests to introduce a bond probability p'f^ < 1 for the nearest neighbour 
particles to belong to the same droplet: 

K = Prob {^-J<o), (7) 

where p is distributed according to a Maxwell law, that is, the probability 
for a particle to have a momentum p is 

One can see that the choice of the bond probability (|^) allows to extend the 
argument to the case in which the loss of a particle of the cluster means the 
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break of several bonds |ni|. We obtain: 

= Prob (p^ < 2mj) = 
Ait 

= 1 - 
= 1 - 



[27rmA;BT]3/2 
4 



1 

oo 

'2mJ 



Prob (p^ > 2mj) 



-p-^/2mkBT 2 



p dp 



(9) 



Performing a Laguerre expansion, one finds [|T^ 



1 - 0.911e"^/^^'s^ - 0.177 I 1 



J 



(10) 



which shows the similarity between the bond probability definition and 
the Coniglio-Klein distribution of probability 



Vh = 1 



-2J/kBT 



We will use this much simpler bond probability to simulate the kinetic en- 
ergy contribution. This definition of bond-correlated percolation has also the 
advantage that the droplets percolate at the thermodynamic transition point 
with the exponents of the Ising model [|1^]. This is not true in general with 
the site-correlated percolation given by the Ising clusters {p^ = 1). 



2.2 Previous results 



The LGM has been introduced in order to describe fragmentation. In refs. 



TT| , |T2[ several analytical simplifications were made on the model, as well 



as numerical simulations taking the analogy with the Ising model at fixed 
magnetic field. The authors used a definition of droplets which lead the 
thermal critical point to coincide with a percolation point. Recently the 
authors of refs. |ll|, O] have also considered isospin-dependent interactions. 



and they have determined the distribution of clusters together with other 
observables 



13]. Similar studies can be found in ref. |20 



In ref. |]TD[ it was noted that there exists a full line, not only a single 
point, in the phase diagram of the LGM, for densities larger than 0.5, where 
distributions of fragment sizes which obey a power law distribution are ob- 
tained. On this line, which is very close to the Kertesz line [^, one finds a 
phase transition because the droplets are defined with a physical prescription. 
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Along this transition, there are no discontinuities for any thermodynamic 
quantities (energy, specific heat,...). However, for some nonlocal quantities, 
as clusters, a discontinuity appears. The critical exponents along this line 
(outside of the critical point) should give some clue about the kinship with 
a percolation transition. 

2.3 Present model: Ising with fixed number 
of particles 

In all the numerical calculations mentioned above, the LGM is equivalent to 
the Ising model with a constant magnetic field. This corresponds to a fixed 
chemical potential. Hence the number of particles (and consequently, the 
density) is conserved on the average in the framework of the grand canonical 
formalism. For systems as small as nuclei, it is certainly appropriate to 
consider the model with a fixed number of particles, which corresponds to 
the canonical ensemble. As it was stated above, the nucleus is a closed system 
without heat bath. Under these conditions, the microcanonical ensemble is 
the most appropriate formalism to use, but we will leave the constraint of the 
energy as an open problem and concentrate on the crucial difference which 
exists between the canonical and the grand canonical in this case. The reason 
is that the constraint of fixed number of particles can not always be satisfied 
with the last formalism, in sharp contrast with the situation in the former 
framework as we shall see. 

The LGM in the canonical formalism is just the Ising model with fixed 
magnetization (IMFM), here with fixed number of particles. We introduce 
the following partition function 

= ^e'^^<'^>'"'''^ ■5(M-M), (12) 

where (3 = J/k^T (we take J/Zcb = 1 for simplicity), M G [—1,1] is a 
fixed number, and the delta function constrains the configurations [a] which 
contribute to the sum to those with fixed density, so that M = M, where M 
is defined in (^) above. 

As already stated above, the grand canonical formalism of the LGM is the 
Ising model with constant h. Of course an Ising model with a magnetic field 
h which produces a magnetization (M) at a certain value of (3 is equivalent 
in the thermodynamic limit (that is, it gives the same value for the different 
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observables) to the model (|T^) with magnetization M = (M) and the same /3. 
This is so because in the hmit V oo, the fluctuations in M are negligible, 
and the distribution is a Dirac delta centered in (M). 

But, reversely, the description of a model with fixed magnetization is 
not necessarily equivalent to a model with fixed magnetic field in the grand 
canonical ensemble: in the LGM (we will refer from now on to the LGM as 
the model considered in the grand canonical formalism), we obtain a certain 
mean magnetization (or density, see eq. (^) with a magnetic field h at some 
value of l3. When we increase /3, we have to lower h in order to get the same 
value of (M). It is always possible to find a value of h which gives a fixed 
(M) while P < jdcih = 0). However, for (3 > Pc, there is an spontaneous 
magnetization {M{f3)), and the presence of a magnetic field h ^ can only 
increase |(M(/?))|, therefore values of |(M)| < |(M(/3))|, are not accessible 
in the grand canonical formalism. In fact, these non-equilibrium states are 
separated from the normal states by the first order transition in /3 mentioned 
in section 2.1. 

However, we can still consider the values of /? and the magnetization that 
correspond to non-equilibrium states in the LGM, in the framework of the 
model (p!2D. They will correspond to a definite mean value of the energy, for 
example. These points of the phase diagram are points of thermodynamic 
equilibrium of the model (0). Hence one may ask about the behaviour of 
the energy when crossing the corresponding transition point of the LGM and 
the existence of a phase transition at this point in the present model. We 
aim to investigate this point in the sequel, as well as its consequence on the 
behaviour of the fragment size distributions. 



2.4 Update algorithm 

The Ising model with constant order parameter (magnetization M) has been 
considered previously from the point of view of its critical dynamic exponents 
(see, eg ||2^), but it seems that its equilibrium properties have not been 



addressed |2^. The corresponding dynamics would be in principle that of 
the Kawasaki spin exchange model [Q, but, as the authors of noted, one 



should take special care with the detailed balance in a numerical algorithm 
for this model. 

In our Monte Carlo computations, the following update algorithm has 
been used: two spins are picked up at random, and they are exchanged 
according to the standard Metropolis prescription 12^] only if they have dif- 
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ferent signs. This algorithm preserves the magnetization (that is, the number 
of spins up, which represent the particles), and it can easily be seen that it 
satisfies detailed balance. 

In fact it is striking to notice how one can arrive to non-Boltzmann dis- 
tributions if one does not care very much about the update procedure. For 
example, we also considered the following algorithm: take one spin sequen- 
tially on the lattice and another one at random, and exchange them only 
if they have different signs. This does not satisfy detailed balance, as the 
following argument shows. Consider an initial configuration with every spin 
up except for a single spin down. Suppose that the sequential process is such 
that we have to take the spin at the position i = 0, which is up. With the 
algorithm described, there is only a possible final configuration, the one in 
which the spin down has been exchanged with the spin up at i = 0. If we 
consider now the reverse process, we have to exchange the spin down at z = 
with a spin up. But now we have many final configurations, not only the ini- 
tial one, so it is clear that this process does not satisfy detailed balance. We 
can modify this sequential algorithm to solve the problem. If the spins taken 
sequentially and randomly are always exchanged, independently of their sign, 
detailed balance is restored. We verified this numerically, so that the results 
with this last algorithm were the same as with the 'full randomly' algorithm 
described above. In fact, an interesting property of this 'random' algorithm 
in which the two spins are taken at random is that here one can choose ex- 
changing the spins in every case, or only when having different signs. Both 
situations verify detailed balance, as we checked numerically, thanks to the 
randomness in the choice of the spins. 

Similar subtleties happen in algorithms where nearest neighbours are ex- 
changed, as in Kawasaki dynamics, which were first noted in p3|. 



3 The (p, T) phase diagram and 
cluster formation 

3.1 Definition of observables 



We simulate the IMFM, defined by eq. ([T^), in a three-dimensional cube of 
V = sites, and with N = p-V spins up and periodic boundary conditions. 
We work out the phase diagram in the {p, P) plane. 
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The thermodynamic observables are the energy 

and its P derivative, the specific heat 

C, = 3V (14) 
The droplet observable which we consider is the second momentum 

S2 = J2s'P{s), (15) 

s 

where P{s) is the number of droplets of size s, and the sum extends over all 
droplets except for the largest one. 

The droplets are defined from the Ising clusters with the probability law 
([TT|). With this definition, the transition properties of the thermodynamic 
transition, given by the divergence of Cy, and the droplet transition, given by 
the divergence of 5*2, coincide in the standard Ising model (constant magnetic 
field) at /i = 0. We want to see whether this is also the case in the IMFM. 

3.2 Phase diagram and cluster formation 

The phase diagram of the IMFM should be the same as the one of the LGM 
because the thermodynamic relations between the different quantities must 
be the same in the domain where the LGM displays an homogeneous phase. 
The thermodynamic transition lines must coincide in both phase diagrams. 
However, the nature of this transition can be different, as noticed above. 

We can also understand physically the existence of a phase transition 
in the Ising model with constant magnetization with the help of Fig. [l|. It 
shows the curve M(T) for the usual Ising model ai h = magnetic field. 
Constraining the system to move in T with fixed magnetization M (indicated 
by the dashed line) will be equivalent to pass through Ising models with 
different values of h (those which give the magnetization M{T) = M at 
each T), until one finally reaches the M{T) curve corresponding to h = 0. 
Then it is not possible to go to the left in the usual Ising model with any 
magnetic field, but this is possible with the present model. One expects a 
phase transition at that point, and this will be confirmed below. Therefore, 
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M 




Figure 1: Plot of M{T) in the Ising model at h = magnetic field, and 
trajectory of constant magnetization. 

the curve M{T) plotted in Fig. |I| will exactly correspond to the transition 
line of this model in the (p, (3) plane. 

In Fig. 1^ we plot the phase diagram obtained directly in a numerical 
simulation of an L = 20 lattice. We represent the thermodynamic transition 
(maximum of the specific heat, divergence in the infinite volume limit) and 
the droplet transition (maximum of S'2). For L = 20 these transitions are 
different, but a priori this could be a finite size effect. We shall see in the next 
sections whether these two transitions are the same. The thermodynamic 
transition is symmetric with respect to p = 0.5, because the thermodynamic 
properties are the same in terms of spins up or down, that is, the model has 
a symmetry a ^ —a which corresponds to p ^ 1 ^ P- The droplet transition 
is not symmetric because we define only droplets for sites corresponding to 
(T = +1. In Fig. 1^ the thermodynamic transition for L = 20 presents an 
"anomaly" at p = 0.5 so that the curve seems not to be concave there. This 
is a finite size effect, due to the fact that the critical exponents vary along 



11 




0.0 0.2 0.4 0.6 0.8 1.0 



P 

Figure 2: Phase diagram of the Ising model with constant magnetization, 
obtained numerically from an L = 20 3(i cubic lattice. The solid line is 
the thermodynamic transition and the dashed line the droplet transition. 
The dotted line indicates the infinite volume limit of the thermodynamic 
transition (see text). The black point at p = 1 is the three-dimensional 
bond-percolation point. The line with the arrow (dot-dashed) symbolizes a 
typical trajectory, see text. 

the line, as we will see in section 4. Then the thermodynamic limit is reached 
differently at different points on the line. Actually, we have included in the 
Figure the infinite volume curve obtained from the values of the critical points 
at p = 0.3 and p = 0.5 calculated in section 4.3. 

The region A in Fig. ^ is the disordered region, B is a region where 
an infinite cluster is formed but the spin system, from a thermodynamic 
point of view, is disordered, and C is a region with an infinite cluster which 
also corresponds to an ordered thermodynamic system. The thermodynamic 
transition is the line which separates the regions B and C, and the droplet 
transition is the one which separates the regions A and B. 

The physical origin of the infinite cluster is, on the one hand, the Ising- 
type interaction, which, in the nuclear context, simulates the strong inter- 
action between the particles. The most probable configuration for large (3 
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is that with all spins up together, forming a spherical cluster. This is only 
relevant in region C. On the other hand, the density also influences the de- 
velopment of clusters, because we know that, even in absence of interaction, 
the system percolates when there are many particles. This effect appears in 
region B, and this is the reason why the droplet transition separates from the 
thermodynamic one for large values of p. Besides that, (3 also influences the 
formation of physical droplets through the bond probability (|llD, simulating 
a kinetic energy for the particles depending on the temperature: the lower 
the temperature, the less the particles can move. 

Far inside the ordered phase, the system favours the configuration which 
minimizes the energy, and then a large cluster is formed. All occupied sites 
inside the cluster contribute positively to the energy, as well as those outside 
the cluster {(JiCJj = +1). The only contributions with (Xjcrj = —1 are those 
located at the border, so that the configuration will show a minimal surface, 
something very similar to a sphere. One can estimate the energy in this case. 
Let us consider a cubic shape for this cluster. Then, the total amount of 
(potential) energy is 

Et = 3F (l - , (16) 

so that 

E^^^l-'J^, (17) 

3V L ' ^ ^ 

and in the thermodynamic limit (L oo), E is normalized to 1. This is why 
we chose the factor 3V in the normalization of E and in eqs. (|r^) and 

dl)- 



3.3 Cluster distribution and nuclear fragmentation 

Along the droplet transition, the observable 5*2 diverges, because there the 
fragment size distribution P{s) follows a power law 

P{s) ~ s~\ (18) 

A possible fragmentation scenario can be suggested in the (p, (3) diagram 
of Fig. 1^ (dot-dashed line). Consider a nuclear collision leading to an ex- 
cited system of particles in a volume V, in thermodynamic equilibrium 
or at least close to it. This point would correspond to a rather large density 
(p ~ 1) and small (3 (large temperature), lying below the droplet transition 
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line. The trajectory would then move towards decreasing p and increasing (3. 
It would cross the critical droplet line at some point (po, /^o) and pursue its 
way towards lower density and temperature, up to some freeze-out where the 
fragment size distribution gets fixed. From thereon the system would con- 
tinue to cool down and expand continuously, crossing the thermodynamic 
transition line at some other point. In practice, of course, the system follows 
a whole set of trajectories corresponding to different experimental initial 
conditions (system, energy, impact parameter,...). Unfortunately it is not 
possible to get direct information about these trajectories. There could exist 
indications that the crossing of the critical droplet line may occur for rather 
high values of the density. In fact, indirect experimental information (agree- 
ment with percolation results) indicates that the exponent r could lie close 
to the value 2.2 which corresponds to p > 0.5 (see section 4.4 below). On the 
other hand, molecular dynamics (MD) calculations seem to indicate that the 
system breaks up rather early during the collision process . Whether the 
critical r can be put in direct correspondence with experimental measure- 
ments of fragment distributions which correspond to averages over different 
events remains however an unsettled question. 



4 Phase transitions and their characteristics 
4.1 Critical exponents 

As we have seen, the phase diagram of the IMFM shows, in principle, two 
phase transitions. We may ask whether these two transitions coincide for 
small values of p, the type of thermodynamic transition and what is its 
order, its universality class, and the value of the exponent r, whether it is 
constant or not along the transition line. 

The critical exponents will bring an answer to these questions. We have 
two observables which diverge at their corresponding transitions in the ther- 
modynamic limit, the specific heat 

C.~(/5-/?J)-", (19) 

which defines the exponent a, and the second momentum 

^2~(/3-/3e°)-^ (20) 

which defines the exponent 7, /3j and (3^ correspond to the thermodynamic 
and droplet transition respectively. Note that in principle, the exponents 
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of the thermodynamic and droplet transitions will be independent, because 
we have two different models: Ising (with fixed magnetization), and a bond- 
percolation correlated model, constructed on the IMFM system. We know 



that in the case of the standard Ising model, the definition (|TTD of the bond 
probability is such that the thermodynamic and the droplet exponents are the 
same, and the position of the transitions coincide. But this could be different 
in the present case. Actually, the exponent 7, defined as the one which gives 
the divergence of the magnetic susceptibility, is not defined for the IMFM, 
because we keep the magnetization fixed as an external parameter. This is 
also the case of the remaining "magnetic" exponents, (3 and rj. Therefore, 
the universality class of the IMFM will be completely determined by a and 
z/, which should obey the hyperscaling relation 

a = 2-ud. (21) 

We will denote the exponents corresponding to the droplet transition with a 
little bar above the corresponding greek letter. 

We can measure critical exponents from standard finite size scaling 
p8| . For finite lattice size L, C„ and S2 show peaks (at different points, Pj{L) 
and Pf{L), respectively), which scale with L as 

C,{L) = A + BL'^'^ (22) 
S2{L) = AV<I\ (23) 

The exponent v can also be obtained from a fit to the law 

/3,(L) - /3,(oo) = (24) 

where /5c(oo) is the thermodynamic value of the transition point, and (3c{L) 
is the value of (3 where the thermodynamically divergent observable has a 
maximum at finite L. 

The exponents of the droplet transition will also obey the corresponding 
scaling and hyperscaling relations 

1 = (25) 

a + 2^ + 7 = 2, (26) 
a = 2-ud, (27) 

where these exponents are defined by the behaviour of various moments of 
the cluster size distributions, see eg [^. 
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Using scaling and hyperscaling with d = 3, the exponent r of the distri- 
bution of fragments can be obtained directly from f/: 

r = 2 + ^ = 2 + l±^ = 2.2 + ^V + oif), (28) 
p + 7 — f] 25 

so that if fj is small (which is usually the case), r is 2.2 modified by a small 
amount. We remark that for all the models with small fj (and for all the 
standard models in three dimensions, f] < 0.1), r can only vary between 2.1 
and 2.3. The only nontrivial result is the deviation from 2.2. 

We shall be able to check the validity of these relations through the deter- 
mination of T by means of two independent methods, directly, by measuring 
the droplet distribution followed by a fit to eq. (plSf), and by computing the 
critical exponents with the use of eq. (^). 

4.2 The numerical simulation 

In order to compute the critical exponents, and to see if the droplet and 
thermodynamic transitions coincide in the infinite volume limit, we have 
studied the transitions with large statistics at three different values of p 
in the phase diagram, p = 0.3,0.5,0.7. The peak in (eq. (P^)) can 
be easily located with the spectral density method ||30|, but unfortunately. 



this is not possible with S2, because we construct this observable with bond 
probabilities for each configuration. This means that the value of S2 is not 
completely determined by the configuration, and then it is not possible to use 
the spectral density method to extrapolate to different values of f3. Therefore, 
we have obtained the maximum of S2 with less precision than in the case of 

Cv 

We studied the thermodynamic transition at p = 0.3 and p = 0.5 (because 
of the symmetry of this transition, p = 0.7 is equivalent to p = 0.3), in order 
to fix the exponents a and u, which inform us about the universality class of 
the IMFM model. In Table |I| we show the lattice sizes and statistics used. 



together with the integrated autocorrelation time |^ for the energy. 

For the droplet transition we only studied lattices up to L = 24, first 
because of the difficulty of finding the peak of 5*2 as mentioned above, and 
second because the exponent r may be a meaningful physical quantity in 
the case of a finite system, which may eventually be directly or indirectly 
confronted with the experiment. 
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p = 0.3 


p = 0.5 


L 


riE) 


Sweeps [xlOV(E)] 


r{E) 


Sweeps [x10^t{E)] 


10 


8.55(3) 


1091 


12.15(6) 


985 


16 


12.84(11) 


932 


22.1(3) 


407 


20 


15.72(12) 


381 


27.5(3) 


345 


24 


20.84(12) 


1061 


35.2(2) 


680 


28 


27.2(2) 


527 


43.7(4) 


738 


32 


33.0(3) 


362 


53.5(7) 


558 


40 


63.9(1.1) 


110 


75.3(8) 


189 


48 


109(2) 


108 


110(7) 


43 



Table 1: Total number of sweeps or update steps (measurements are per- 
formed every 10-15 sweeps) and the corresponding integrated correlation 
times r for E, used in the analysis of the thermodynamic transition at p = 0.3 
and p = 0.5. 



In the simulation we proceed as follows: we make 10-15 update steps (in 
each step the algorithm described in section 2.4 is applied V times), and 
then we perform a measurement. At this moment the droplets are defined 
according to (^)- Hence we can calculate 82- The errors are computed with 



the jackknife method 32 



4.3 The thermodynamic transition 

We start by considering the nature of the thermodynamic phase transition 
(solid line in Fig. 0), ie the points in the phase diagram where the specific 
heat diverges. We show in Fig. | the caloric curve {E vs. (3) for two lattice 
sizes, L = 10, 24. In principle, these kind of curves are not very precise 
to distinguish the order of the phase transition. A first order transition 
means a discontinuity in this curve, and a second order one means a vertical 
slope at the critical line. But this is only true in the thermodynamic limit, 
L = 00. In a finite lattice there are not singularities at all, and then it 
is necessary to carry out a careful study of the thermodynamic limit by 
considering finite size scaling techniques, such as expressions ( ^21) and (^31). 
The values of the critical exponents obtained in this way inform us about the 
order of the phase transition, and define the universality class in the case of 
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Figure 3: Caloric curve {E vs. (3) for two finite lattice sizes, L = 10 and 
L = 24, at p = 0.3. 

a second order transition. We will show evidences of the second order nature 
of the transition line at every value of p. On the one hand, a single peak 
in the histogram of the energy was observed in all numerical simulations. 
This is a strong indication of a second order phase transition. On the other 
hand, this was confirmed by an analysis of the critical exponents u and a, 
which show a different behaviour from that of a first order transition (a = 1, 

iy=l/d= 1/3). 

Exponents u and a determine the universality class of the transition. In 
order to see the values of these exponents, we have to find the behaviour in 
the V ^ oo limit. We consider lattice sizes up to L = 48 and study these 
exponents by finite size scaling at p = 0.3 and p = 0.5. 

In Table ^ we report the values of the thermodynamic transitions and 
the maxima of for the different sizes at p = 0.3 and p = 0.5. 

Consider first the case p = 0.3. One needs to make a fit to eq. (p2D 
in order to obtain a/u. If the constant term is not very important, a plot 
InC™^^ vs. InL should fit to a straight line of slope a/u. This plot is shown 
in Fig. ^ where we see that the slope changes with L for small L. This means 
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p = 0.3 


p = 0.5 


L 








^max 


10 


0.2698(3) 


3.063(5) 


0.2778(10) 


3.880(7) 


16 


0.2516(3) 


5.109(9) 


0.2495(4) 


6.380(14) 


20 


0.2445(2) 


6.191(8) 


0.2417(2) 


7.739(11) 


24 


0.23988(11) 


7.033(7) 


0.23694(6) 


8.819(9) 


28 


0.23658(11) 


7.758(8) 


0.23364(7) 


9.726(11) 


32 


0.23435(11) 


8.45(2) 


0.23157(8) 


10.519(12) 


40 


0.23152(5) 


9.72(2) 


0.22862(7) 


11.80(3) 


48 


0.229885(14) 


11.11(4) 


0.22680(8) 


12.93(6) 



Table 2: Thermodynamic transition /?J and the maxima of for different 
lattice sizes at p = 0.3 and p = 0.5. 

that corrections to scaling are important. We observe that the slope reaches 
an asymptotic limit from L = 24 on. Actually, the fit for L = 24,28,32,40 
is very good: 

a/u = 0.634(4) xV'DF = 0.290/2, (29) 

where DF stands for the number of degrees of freedom in the fit. The L = 48 
point, however, stays a bit away of this value. In fact, if we try to include 
this point, the fit obtained is not good because we obtain a large x^/DF, 
and it gives a value of a/u = 0.643(4). It is possible that the failure of this 
fit is due to a subestimate of the error in this last point, which has a large 
autocorrelation time (see Table |1]). To be cautious, anyway, we will take as 
the best estimate for a/u, the value 

a/u = 0.639(8). (30) 

An independent fit from this one is the two-parameter fit to eq. (^4]), 
which gives the best fit for u = 0.60(4). Here the error is very large. Consid- 
ering instead the result (pOD, hyperscaling (pT]) gives for u and a the values 

u = 0.5496(12) a = 0.351(4). (31) 

One can check this value of u by means of a fit to eq. (^4]) . Taking again the 
lattices 24 < L < 48, we obtain 

/?c(oo, p = 0.3) = 0.22606(13) xVDF = 0.881/2, (32) 



19 



2.5 — 



2.0 



1.5 



p = 0.3 
p = 0.5 
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Figure 4: Plot of InC^^^ vs. InL for L = 10,16,20,24,28,32,40,48 at 
p = 0.3 and p = 0.5. A straight line fit should give a/u. We observe scaling 
corrections for small L in both cases. 

which is a good fit. If we add the L = 48 lattice, this value is slightly modified 
iPc = 0.22590(6)), and the xVDF is still acceptable (3.5/3). Therefore, we 
will take the values ( |3TD as the best estimate for a and u at p = 0.3. 

Next, consider the case p = 0.5. Here a plot \nC^^^ vs. InL does not 
give a good fit up to L = 40, see Fig. ^, but if we include the simulation for 
L = 48, it seems that an asymptotic behaviour is reached for the three last 
lattices L = 32,40,48. The fit gives 

a/z/ = 0.513(10) xVDF = 0.05/1, (33) 

which, using hyperscaling, means 

u = 0.569(2) a = 0.292(6). (34) 

A fit to eq. (|2^) (using the three last lattice sizes) with this value of z/ is not 
bad. It gives 

/5c(oo, p = 0.5) = 0.2223(3) xVDF = 0.895/1, (35) 
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However, this value of /?c(oo) is not compatible with the one we had expected, 
the Ising value (3c{oo) ~ 0.22165, though it is not far from it. In fact, if we 
consider the Ising value of exponent u ~ 0.63, and use it to obtain /3c(oo) 
in the fit to eq. (0), we obtain quite a good fit with a value of Pc{oo) 
compatible with the Ising one: 

Pcioo) = 0.2216(2) xVDF = 0.274/1, (36) 

This means that the behaviour observed for L = 32, 40, 48 can be still transi- 
tory and we should not discard the three-dimensional Ising values at p = 0.5. 

The exponents which we determined were obtained with an accuracy 
which is high enough to conclude that their variation along the transition 
line is soft. This requested a substantial numerical effort in the framework 
of the canonical ensemble. Working in the framework of the microcanonical 
ensemble may request less numerical investments. A more detailed study, 
using different values of p, should be done in order to clarify the universality 
classes along the transition line. 

4.4 The droplet transition 

We consider now the droplet transition, ie the points of the phase diagram 
where 5*2 diverges, which means that a droplet percolates through the system. 
For p < 0.5 this transition line is very close to the thermodynamic one, but for 
p > 0.5, it seems to be very different, and strictly under the thermodynamic 
line, which is easy to understand because it is the high density of particles 
and not the interaction which drives the system to percolation in this case. 
One expects that for p < 0.5, the two transition lines coincide. We can check 
this, as well as the value of exponent r at different points of the droplet line. 

In Table ^, we report the l3^{L) values of the droplet transition and the 
maxima of 5*2 for different lattice sizes atp = 0.3, p = 0.5 and p = 0.7. 

A fit of the maxima of 5*2 to eq. ( [23| ) gives the exponent 7/z/, and there- 
fore, f] from eq. (|25|), which will allow us to obtain r through eq. (pS]) . 

We use only data from lattices of sizes L = 10, 16, 20, 24, see Table |^, 
because of the difficulty of locating the droplet transition (as mentioned 
before, it is not possible to use the spectral density method in this case). 
With these sizes, we are not able to obtain a stable value of j/u which can 
be extrapolated to the infinite volume. We will simply obtain estimates for 
the exponent r which will depend on the lattice size. These are shown in 
Table |. 
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p = 0.3 


p = 0.5 


p = 0.7 


L 








^max 






10 


0.2560(5) 


5.149(2) 


0.225(3) 


8.62(2) 


0.194(2) 


11.70(5) 


16 


0.2440(2) 


9.402(12) 


0.2230(5) 


20.36(4) 


0.1984(2) 


30.43(5) 


20 


0.23960(10) 


11.91(2) 


0.2227(4) 


30.58(4) 


0.1990(6) 


48.15(5) 


24 


0.2367(3) 


14.16(2) 


0.2227(2) 


42.98(4) 


0.2005(6) 


70.6(2) 



Table 3: Droplet transition (3^ and the maxima of 5*2 for different lattice 
sizes at p = 0.3, p = 0.5 and p = 0.7. 





L = 10,16 


L = 16,20 


L = 20,24 


p 




r 


7/z/ 


r 


7/z/ 


r 


0.3 


1.281(3) 


2.4015(12) 


1.060(9) 


2.478(3) 


0.949(12) 


2.519(5) 


0.5 


1.829(6) 


2.242(3) 


1.823(11) 


2.244(3) 


1.867(9) 


2.233(2) 


0.7 


2.034(10) 


2.192(2) 


2.056(9) 


2.188(2) 


2.10(2) 


2.176(5) 



Table 4: Values of 7/z/ and r, taking the data from Table We observe the 
evolution with p for every lattice size. 



We can check now the values of r with the second method: directly from 
the histogram of the number of droplets of a certain size. In Fig. |^ we show 
some of these histograms, for L = 20 at p = 0.3, p = 0.5 and p = 0.7, in 
logarithmic scale. The actual value of the slopes are not very precise, because 
they depend on the region where the fit is performed. Nevertheless, there 
are central zones in every histogram which give slopes stable under small 
variations of their length, and where we retrieve the values of r obtained in 
Table ^. This means that scaling and hyperscaling are satisfied, and we see 
that the results show a dependence of exponent r on p (besides the finite size 
dependence, of course), which suggests a continuous change of the values of 
this exponent along the droplet transition line. 

Finally, one may raise the question of whether the transitions merge to- 
gether at the thermodynamic limit. We show the differences between (3j and 
/?° at p = 0.3, 0.5, 0.7 for the lattice sizes L = 16, 20, 24 in Table |. 

If we try a fit of the data of Table ^ to Af3 = AL~^, we obtain: 

p = 0.3 a; = 2.1(2) = 0.48(5), xV^^ = 0-28/1 (37) 
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Figure 5: Number of droplets P{s) with s particles in logarithmic scale. The 
slopes of the fitted straight lines give exponent r. This is done for L = 20 at 
p = 0.3,0.5,0.7. 



p = 0.5 x = 1.54(6) x~^ = 0.65(3), x^/DF = 0.19/1 (38) 
p = 0.7 x = 0.73(4) x-^ = 1.37(7), x^/DF = 0.31/1 (39) 

We see that the fits (|37D and (|38| ) for p = 0.3 and p = 0.5 respectively, 
are quite good. One could argue that if exponents u and z/ were the same, 
this value should be equal to in the fit above. In fact we see that in (^Tf), 
the value of is near to (although not completely compatible with) the 
value of u that we obtained in eq. (^). In the case p = 0.5, we obtain a 
value of x~^ which is compatible with the 3(i-Ising u exponent (0.63). This 
agrees with our conclusion at the end of the last section that a compatibility 
with the Ising exponents should not be discarded at p = 0.5. 

Finally, we see from (|39|) that a fit to A/3 = in the L ^ oo limit is also 
possible with these data. However, the value of x is now less than 1 and has 
no relation with exponent u. We know that the droplet transition line must 
end at the bond-percolation point (see Fig. |^), so that we expect that this fit 
is a spurious effect, or, equivalently, that exponent x in this kind of fit will 
approach zero for larger lattice sizes. We then conclude that for p > 0.5 the 
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L 


A/?(0.3) 


A/3(0.5) 


A/3(0.7) 


16 


0.0076(4) 


0.0265(6) 


0.0532(4) 


20 


0.0049(2) 


0.0190(4) 


0.0455(6) 


24 


0.0032(3) 


0.0142(2) 


0.0394(6) 



Tabic 5: Differences between the thermodynamic and the droplet transition 
points at p = 0.3, 0.5, 0.7. 

thermodynamic and droplet transitions are different. 

5 Conclusions 

In the present work we presented and analyzed a simple though realistic 
2>d Ising-type model with fixed magnetization (density) which we considered 
as a generic description of nuclear matter fragmentation in thermodynamic 
equilibrium. We worked out the phase diagram of the system in terms of 
density and temperature. We analyzed the thermodynamic properties of the 
system as well as the properties of the size distribution of bound fragments. 
For compact systems (relative density larger than 0.5) we found two types 
of transitions, one corresponding to a percolation transition for the system 
which breaks into pieces, and a thermodynamic transition at lower temper- 
ature. For dilute systems (relative density lower than 0.3), which should 
correspond to situations where the expanding system reaches the freeze-out 
point, the two transitions merge and the exponents related to the thermody- 
namics and the fragment distributions change with the density. 

There are strong indications that the scaling and hyperscaling relations 
between the exponents are satisfied. In all cases which were worked out 
numerical results indicate that the transitions can be interpreted as being 
second order. For this reason, one cannot expect to observe the experimen- 
tally controversed latent heat plateau in the caloric curve, at least not in 
a clear-cut way in the small systems which we considered here. One may 
perhaps argue here that the energy fiuctuations which are present in the 
canonical treatment performed here smear out the energy and hence pre- 
clude a clear-cut observation of the plateau. We insist however on the fact 
that the critical exponents definitely tell us that the transitions are second 
order as already mentioned above even if at the present stage of our investiga- 
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tions we cannot say with certitude to which universahty class the transitions 
belong (it seems clear, however, that they experiment a soft variation along 
the transition line). Actually, the fact that the transition is first order for 
non-fixed magnetization, and second order when the magnetization is fixed 
could have been expected somehow. For a fixed magnetization, the configu- 
ration space is smaller, there are less fluctuations, and it is natural that the 
transition gets weaker. 

It remains nevertheless true that the fragmenting nuclear system is iso- 
lated (if one considers the whole system which enters the reaction process 
leading to its fragmentation) which justifies a microcanonical treatment as 
advertised by D. Gross and collaborators [§. It may be worthwhile to con- 
sider this point, although it is not clear to us whether it would clarify the 
situation since in the microcanonical ensemble it is the temperature which 
fluctuates strongly at the transition points for fixed energy, when the reverse 
is true in the canonical ensemble. It is our aim to work out the present model 
in this framework in the next future. 

Last but not least, the present work shows clearly that models which 
describe finite size systems and which differ by seemingly harmless charac- 
teristics (strictly fixed number of particles in IMFM here, fixed number of 
particles in the average in the LGM) can lead to quite different results as 
far as the order of the transition is concerned. Experimental data should 
be enough precise and complete in order to be able to distinguish between 
details like the order of the transition. This is certainly a great challenge and 
an incentive to improve our information about the properties of fragmenting 
nuclei. 
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